home *** CD-ROM | disk | FTP | other *** search
/ Meeting Pearls 1 / Meeting Pearls Vol 1 (1994).iso / installed_progs / gfx / lise2.1 / lise / src / b_hf.c < prev    next >
Encoding:
C/C++ Source or Header  |  1993-03-31  |  15.1 KB  |  685 lines
  1. /* standard begin
  2. */
  3.  
  4. #define MAXTERM 30
  5.  
  6. #include <stdio.h>
  7. #ifdef AMIGA
  8. #include <gfxamiga.h>
  9. #include <libraries/dos.h>
  10. #include <libraries/dosextens.h>
  11. #else
  12. #include <gfx.h>
  13. #endif
  14.  
  15. float  *spc,*err,*tim;
  16.  
  17. double x,y;
  18. double Pi = 3.1415927,
  19.        Gs = 2.0,
  20.        k = 1.38062E-23,
  21.        Mu0 = 1.256637E-6,
  22.        MuB = 9.2741E-24,
  23.        gI = 0.06,
  24.        Nq = 4.0,
  25.        l = 3.0,
  26.        Nel = 5.0,
  27.        J_grund = 0.0,
  28.        Fak_BN = 0.0,  /* must be set within programm */
  29.        R_mittel = 0.0, /* must be set within programm */
  30.        B_core = 10.0,
  31.        Tbeg = 50.0,
  32.        Tend = 1000.0,
  33.        Tstep = 10.0;
  34.  
  35. double *EJf,
  36.        *Lges,
  37.        *Sges,
  38.        *Jges;
  39.  
  40. double *sum1,
  41.        *sum2,
  42.        *sum3;
  43.  
  44. double CurieTerm,
  45.        VanVleckTerm;
  46.  
  47. int    NTerm = 0,
  48.        Flag_1dT = 1;
  49.  
  50. char   *fnam;
  51.  
  52. help()
  53. {
  54.   printf("this programm calculates the paramagnetic enhancement factor for\n");
  55.   printf("various ions at various temperatures.\n");
  56.   printf("B_HF [-r file] [-T] [-Ts n.m] [-Te n.m] [-Ti n.m] [-go]\n");
  57.   printf("      -r file    reads data from file\n");
  58.   printf("      -T         generates spectrum as function of direct Temperature\n");
  59.   printf("      -Ts n.m    sets start temperature\n");
  60.   printf("      -Te n.m    sets end temperature\n");
  61.   printf("      -Ti n.m    sets temperature increment\n");
  62.   printf("      -go        suppress interactive mode\n");
  63. }
  64.  
  65. double dabs(x)
  66. double x;
  67. {
  68.    if(x >= 0.0) return(x);
  69.    return(-1.0 * x);
  70. }
  71.  
  72. double Energy(E,us,dir)
  73. double E;
  74. char *us;
  75. int dir;
  76. {
  77. double c;
  78.  
  79.    c = 0.0;
  80.    if(strcmp(us,"J") == 0) c = 1.0;
  81.    if(strcmp(us,"K") == 0) c = 1.38062E-23;
  82.    if(strcmp(us,"eV") == 0) c = 1.60219E-19;
  83.    if(strcmp(us,"erg") == 0) c = 1.0E-7;
  84.    if(c == 0.0) {
  85.       fprintf(stderr,"unrecognised Energy unit :%s\n",us);
  86.       return(E);
  87.    }
  88.    if(dir < 0) c = 1/c;
  89.    return(c * E);
  90. }
  91.  
  92. double ReCuR(R,us,dir)
  93. double R;
  94. char *us;
  95. int dir;
  96. {
  97. double c;
  98.  
  99.    c = 0.0;
  100.    if(strcmp(us,"m") == 0) c = 1.0;
  101.    if(strcmp(us,"cm") == 0) c = 1.0E6;
  102.    if(strcmp(us,"A") == 0) c = 1.0E30;
  103.    if(strcmp(us,"a.u.") == 0) c = 6.748321E30;
  104.    if(c == 0.0) {
  105.       fprintf(stderr,"unrecognised radius unit :%s\n",us);
  106.       return(R);
  107.    }
  108.    if(dir < 0) c = 1/c;
  109.    return(c * R);
  110. }
  111.  
  112. double Tesla(E,us,dir)
  113. double E;
  114. char *us;
  115. int dir;
  116. {
  117. double c;
  118.  
  119.    c = 0.0;
  120.    if(strcmp(us,"T") == 0) c = 1.0;
  121.    if(strcmp(us,"KG") == 0) c = 0.1;
  122.    if(c == 0.0) {
  123.       fprintf(stderr,"unrecognised Magnetic Field unit :%s\n",us);
  124.       return(E);
  125.    }
  126.    if(dir < 0) c = 1/c;
  127.    return(c * E);
  128. }
  129.  
  130.  
  131. double v(l,Lges,Sges)
  132. double l,Lges,Sges;
  133. {
  134. double erg;
  135.  
  136.    erg = (4.0 * Sges - (2.0 * l + 1)) / ((2.0 * l - 1.0) * ( 2.0 * l + 3.0) * Sges * ( 2.0 * Lges - 1.0));
  137.    return(erg);
  138. }
  139.  
  140. double meJNJ(l,Lges,Sges,J)
  141. double l,Lges,Sges,J;
  142. {
  143. double erg,JJ,LL,SS;
  144.  
  145.    JJ = J * (J + 1.0);
  146.    LL = Lges * (Lges + 1.0);
  147.    SS = Sges * (Sges + 1.0);
  148.    if(J == 0) {
  149.       erg = 1.0 + v(l,Lges,Sges) * 0.5 * Gs * (3.0 * JJ - 2.0 * LL - 3.0 * SS - 1.0);
  150.    } else {
  151.       erg = 3.0 * JJ - 2.0 * LL - 3.0 * SS - (3.0 * JJ - LL + SS) * (3.0 * JJ - LL - 3.0 * SS) / (4.0 * JJ);
  152.       erg = erg * v(l,Lges,Sges) * 0.5 * Gs;
  153.       erg = erg + (JJ + LL - SS)/(2.0 * JJ);
  154.    }
  155.    return(erg);
  156. }
  157.  
  158. double sqr(x)
  159. double x;
  160. {
  161.    return(x * x);
  162. }
  163.  
  164. double meJAJ(Lges,Sges,J)
  165. double Lges,Sges,J;
  166. {
  167. double erg,JJ,LL,SS;
  168.  
  169.    JJ = J * (J + 1.0);
  170.    LL = Lges * (Lges + 1.0);
  171.    SS = Sges * (Sges + 1.0);
  172.  
  173.    if(J == 0.0) return(Gs);
  174.    erg = 1.0 + (Gs - 1) * (JJ + SS - LL) / (2.0 * JJ);
  175.    return(erg);
  176. }
  177.  
  178. double Bj(l,Lges,Sges,J)
  179. double l,Lges,Sges,J;
  180. {
  181. double erg;
  182.  
  183.    erg = (sqr(Sges + Lges + 1.0) - J * J) * (J * J - sqr(Lges - Sges))/(6.0 * J);
  184.    erg = erg * (1.0 + v(l,Lges,Sges) * (3 * J * J - Lges * (Lges + 1.0) - 3.0 * Sges * (Sges + 1.0)) / 2.0);
  185.    return(erg);
  186. }
  187.  
  188. double Calc_Bc(gj,J)
  189. double gj,J;
  190. {
  191. return(-1.0 * B_core * (gj - 1) * J);
  192. }
  193.  
  194. calc_ELSJ()
  195. {
  196. double Sg,Lg,x1,x2;
  197. int i;
  198.  
  199.    Sg = 0.5 * ((2.0 * l + 1.0) - dabs(2.0 * l + 1.0 - Nel));
  200.    Lg = Sg * dabs(2.0 * l + 1.0 - Nel);
  201.    J_grund = Sg * dabs(2.0 * l - Nel);
  202.    Fak_BN = 2.0 * MuB * Mu0 * R_mittel / (4.0 * Pi);
  203.    for(i = 0; i < MAXTERM; i++) {
  204.       Lges[i] = Lg;
  205.       Sges[i] = Sg;
  206.    }
  207.    x1 = dabs(Lg - Sg);
  208.    x2 = Lg + Sg;
  209.    NTerm = (int) (x2 - x1 + 1.00001);
  210.    for(i = 0; i < NTerm; i++) Jges[i] = (x1 + (double)i);
  211. }
  212.  
  213.  
  214. double Bi0(J)
  215. double J;
  216. {
  217. double erg;
  218. int n;
  219.  
  220.    n = 0;
  221.  
  222.    erg = Fak_BN * meJNJ(l,Lges[n],Sges[n],J) * J;
  223.    erg = erg + Calc_Bc(meJAJ(Lges[n],Sges[n],J),J);
  224.    return(erg);
  225. }
  226.  
  227. double dexp(x)
  228. double x;
  229. {
  230. double n,fak,px,erg,last;
  231.  
  232.    fak = 1.0; px = 1.0; n = 0.0; erg = 1.0; last = 0.0;
  233.    while(erg != last) {
  234.       last = erg;
  235.       n = n + 1.0; fak = fak * n; px = px * x;
  236.       erg = erg + px / fak;
  237.       if(n > 100.0) break;
  238.    }
  239.    return(erg);
  240. }
  241.  
  242. double Ehoch(x)
  243. double x;
  244. {
  245.    if(x > 0.0) {
  246.       if(x > 500.0) return(dexp(500.0)); else return(dexp(x));
  247.    } else {
  248.       if(x < -500.0) return(0.0); else return(dexp(x));
  249.    }
  250. }
  251.  
  252. /* initialize sum1, sum2, sum3 */
  253. sum_init()
  254. {
  255. double J, gj, JBJ, JBJ1, f1, f2, Sg, Lg;
  256. double dif;
  257. int i;
  258.  
  259.    for(i = 0; i < MAXTERM; i++) {
  260.       sum1[i] = 0.0;
  261.       sum2[i] = 0.0;
  262.       sum3[i] = 0.0;
  263.    }
  264.  
  265.    for(i = 0; i < NTerm; i++) {
  266.  
  267.       J = Jges[i]; Lg = Lges[i]; Sg = Sges[i];
  268.  
  269.       gj = meJAJ(Lg,Sg,J);
  270.       JBJ = Fak_BN * meJNJ(l,Lg,Sg,J);
  271.       f1 = (2.0 * J + 1.0) * gj * MuB * (J + 1.0) * (JBJ * J + Calc_Bc(gj,J)) / (3.0 * k);
  272.       sum1[i]= f1;
  273.  
  274.       if((i > 0) && (J != 0.0)) {
  275.          JBJ1 = Fak_BN * Bj(l,Lg,Sg,J);
  276.          dif = EJf[i] - EJf[i-1];
  277.          f2 = MuB * JBJ1 / dif;
  278.          sum2[i] = f2;
  279.       }
  280.  
  281.       sum3[i] = (2.0 * J  + 1.0);
  282.    }
  283. }
  284.  
  285. double Beta(T)
  286. double T;
  287. {
  288. double s1, s2, s3, kT, EXP_EJdkT;
  289. int i;
  290.  
  291.    kT = k * T; s1 = 0.0; s2 = 0.0; s3 = 0.0;
  292.    for(i = 0; i < NTerm; i++) {
  293.       EXP_EJdkT = Ehoch(-1.0 * EJf[i] / kT);
  294. /*
  295. printf("i = %d , EJ/kt = %f , Exp(EJ/kT) = %f sum1[i] = %f\n",i,EJf[i]/kT,EXP_EJdkT,sum1[i]);
  296. */
  297.       s1 = s1 + sum1[i] * EXP_EJdkT / T;
  298.       if(i > 0) {
  299.          s2 = s2 + sum2[i] * (Ehoch(-1.0 * EJf[i-1] / kT) - EXP_EJdkT);
  300. /*
  301.          printf("s2 = %-E , sum2[%d] = %-E \n",s2,i,sum2[i]);
  302. */
  303.       }
  304.       s3 = s3 + sum3[i] * EXP_EJdkT;
  305.    }
  306. /*
  307. printf("s1 = %f , s2 = %f , s3 = %f\n",s1,s2,s3);
  308. */
  309.    CurieTerm = 1.0E99;
  310.    VanVleckTerm = 1.0E99;
  311.    if(s3 != 0.0) {
  312.       CurieTerm = s1 / s3;
  313.       VanVleckTerm = s2 / s3;
  314.       return(1.0 + (s1 - s2) / s3);
  315.    } else {
  316.       return(1.0E99);
  317.    }
  318. }
  319.  
  320. generate_spectrum(flag)
  321. int flag;
  322. {
  323. int i;
  324. double T;
  325. float *curie, *vanvleck;
  326.  
  327.    if(flag == 2) {
  328.       curie = (float *) calloc(_MAXSPCLEN,sizeof(float));
  329.       vanvleck = (float *) calloc(_MAXSPCLEN,sizeof(float));
  330.    }
  331.    i = 0; T = Tbeg;
  332.    while(T < Tend) {
  333.       spc[i] = (float) Beta(T);
  334.       err[i] = 0.0;
  335.       tim[i] = (float) T; if(Flag_1dT == 1) tim[i] = (float) (1.0 / T);
  336.       if(flag == 2) {
  337.          curie[i] = CurieTerm;
  338.          vanvleck[i] = VanVleckTerm;
  339.       }
  340.       T = T + Tstep; i = i + 1;
  341.       if(i > _MAXSPCLEN) break;
  342.    }
  343.    unlink("Curie.spc");
  344.    unlink("Curie.err");
  345.    unlink("Curie.tim");
  346.    unlink("VanVleck.spc");
  347.    unlink("VanVleck.err");
  348.    unlink("VanVleck.tim");
  349.    writespec("BHF_out",spc,err,i,2,"B_HF-output");
  350.    writespec("BHF_out.tim",tim,err,i,2,"B_HF-output");
  351.    if(flag == 2) {
  352.       writespec("Curie",curie,err,i,2,"Curie Term");
  353.       writespec("Curie.tim",tim,err,i,2,"Curie Term");
  354.       writespec("VanVleck",vanvleck,err,i,2,"Van Vleck Term");
  355.       writespec("VanVleck.tim",tim,err,i,2,"Van Vleck Term");
  356.       free(curie); free(vanvleck);
  357.    }
  358. }
  359.  
  360. plot_spectrum()
  361. {
  362. FILE *fp;
  363.  
  364.    SysCall("lise:ash bhf_out -v -pag");
  365.    fp = fopen("Curie.spc","r");
  366.    if(fp != NULL) {
  367.       SysCall("lise:ash Curie -p 1 -1 -v -pag");
  368.       fclose(fp);
  369.    }
  370.    fp = fopen("VanVleck.spc","r");
  371.    if(fp != NULL) {
  372.       SysCall("lise:ash VanVleck -p 1 -2 -v");
  373.       fclose(fp);
  374.    }
  375. }
  376.  
  377.  
  378. /* interactive Beta table */
  379. table()
  380. {
  381. char *s;
  382. int i;
  383. double T,y;
  384.  
  385.    s = (char *) malloc(100);
  386.    while(1 == 1) {
  387.       printf("next temperature (or Just return to leave):"); fflush(stdout);
  388.       fgets(s,80,stdin);
  389.       if(s[0] == 10) break;
  390.       sscanf(s,"%lf",&T);
  391.       y = Beta(T);
  392.       printf("CurieTerm = %-E , VanVleckTerm = %-E\n",CurieTerm,VanVleckTerm);
  393.       printf("ß(%f) = %-E\n",T,y);
  394.    }
  395.    free(s);
  396. }
  397.  
  398. show_faktor()
  399. {
  400. double J, gj, JBJ, JBJ1, f1, f2, Sg, Lg;
  401. double dif;
  402. int i;
  403.  
  404.    for(i = 0; i < MAXTERM; i++) {
  405.       sum1[i] = 0.0;
  406.       sum2[i] = 0.0;
  407.       sum3[i] = 0.0;
  408.    }
  409.  
  410.    for(i = 0; i < NTerm; i++) {
  411.  
  412.       J = Jges[i]; Lg = Lges[i]; Sg = Sges[i];
  413.  
  414.       gj = meJAJ(Lg,Sg,J);
  415.       printf("gj(L=%2.2f,S=%2.2f,J=%2.2f) = %f\n",Lg,Sg,J,gj);
  416.       JBJ = meJNJ(l,Lg,Sg,J);
  417.       printf("<J|N|J>(l=%2.2f,L=%2.2f,S=%2.2f,J=%2.2f) = %f\n",l,Lg,Sg,J,JBJ);
  418.       JBJ = Fak_BN * JBJ;
  419.       f1 = (2.0 * J + 1.0) * gj * MuB * (J + 1.0) * (JBJ * J + Calc_Bc(gj,J)) / (3.0 * k);
  420.  
  421.       if((i > 0) && (J != 0.0)) {
  422.          JBJ1 = Bj(l,Lg,Sg,J);
  423.          printf("gj-1(l=%2.2f,L=%2.2f,S=%2.2f,J=%2.2f) = %f\n",l,Lg,Sg,J,JBJ1);
  424.          JBJ1 = Fak_BN * JBJ1;
  425.          dif = EJf[i] - EJf[i-1];
  426.          f2 = MuB * JBJ1 / dif;
  427. /*
  428. printf("EJ[%d] = %-E , dif = %-E , JBJ1 = %-E\n",i,EJf[i],dif,JBJ1);
  429. */
  430.       }
  431.    }
  432. }
  433.  
  434. double edit_double(text,x)
  435. char *text;
  436. double x;
  437. {
  438. char *s;
  439. double y;
  440.  
  441.    s = (char *) malloc(100);
  442.    printf("%s = %f ?",text,x); fflush(stdout);
  443.    fgets(s,80,stdin);
  444.    y = x;
  445.    if(s[0] != 10) sscanf(s,"%lf",&y);
  446.    free(s);
  447.    return(y);
  448. }
  449.  
  450. ask_data()
  451. {
  452. int i,n;
  453. double x,y;
  454. char s[4];
  455.  
  456. /*   gI = edit_double("Kern g Faktor",gI); */
  457.    Nq = edit_double("Hauptquantenzahl",Nq);
  458.    l = edit_double("Bahndrehimpuls (s=0,p=1,d=2,f=3)",l);
  459.    Nel = edit_double("Anzahl der Elektronen in der Schale",Nel);
  460.    x = ReCuR(R_mittel,"a.u.",-1);
  461.    x = edit_double("<r-3> in atomaren Einheiten",x);
  462.    R_mittel = ReCuR(x,"a.u.",1);
  463.    x = Tesla(B_core,"T",-1);
  464.    x = edit_double("Corepolarisation in Tesla",x);
  465.    B_core = Tesla(x,"T",1);
  466.  
  467.    printf("calculating defaults ? (Y/N)"); fflush(stdout);
  468.    fgets(s,2,stdin);
  469.    if((s[0] == 'y') || (s[0] == 'Y')) calc_ELSJ();
  470.  
  471.    x = (double)NTerm;
  472.    NTerm = (int) (0.00001 + edit_double("Anzahl der Terme",x));
  473.    for(i = 0; i < NTerm; i++) {
  474.       x = Energy(EJf[i],"K",-1);
  475.       printf("%d : J = %f   S = %f   L = %f   EJ = %f\n",i,Jges[i],Sges[i],Lges[i],x);
  476.       Jges[i] = edit_double("J",Jges[i]);
  477.       Lges[i] = edit_double("L",Lges[i]);
  478.       Sges[i] = edit_double("S",Sges[i]);
  479.       x = edit_double("EJ in Kelvin",x);
  480.       EJf[i] = Energy(x,"K",1);
  481.    }
  482.    printf("pre calculating sums...\n");
  483.    sum_init();
  484.    printf("ready\n");
  485. }
  486.  
  487. calculate_data()
  488. {
  489. int i,n;
  490. double x,y;
  491. int flag;
  492. char c,s[4];
  493.  
  494.    Tbeg = edit_double("Tbeg (Kelvin)",Tbeg);
  495.    Tend = edit_double("Tend",Tend);
  496.    Tstep = edit_double("Tstep",Tstep);
  497.    flag = 1;
  498.    printf("generate Curie curve and Van Vleck curve separately ? (Y/N)"); fflush(stdout);
  499.    fgets(s,2,stdin); c = s[0];
  500.    if((c == 'y') || (c == 'Y')) flag = 2;
  501.    generate_spectrum(flag);
  502. }
  503.  
  504. save_data()
  505. {
  506. FILE *fp;
  507. int i,n;
  508. double x;
  509. char *s;
  510.  
  511.    s = (char *) malloc(100);
  512.    printf("saving to >%s< ?",fnam); fflush(stdout);
  513.    fgets(s,80,stdin);
  514.    if(s[0] != 10) {s[strlen(s)-1] = 0; strcpy(fnam,s);}
  515.    free(s);
  516.    
  517.    fp = fopen(fnam,"w");
  518.    if(fp == NULL) {
  519.       fprintf(stderr,"could not open file for write: %s\n",fnam);
  520.       return();
  521.    }
  522.  
  523.    fprintf(fp,"%f\n",gI);
  524.    fprintf(fp,"%f\n",Nq);
  525.    fprintf(fp,"%f\n",l);
  526.    fprintf(fp,"%f\n",Nel);
  527.    fprintf(fp,"%f\n",J_grund);
  528.    fprintf(fp,"%f\n",R_mittel);
  529.    fprintf(fp,"%f\n",B_core);
  530.    fprintf(fp,"%d\n",NTerm);
  531.    for(i = 0; i < MAXTERM; i++) {
  532.       fprintf(fp,"%f\n",Jges[i]);
  533.       fprintf(fp,"%f\n",Lges[i]);
  534.       fprintf(fp,"%f\n",Sges[i]);
  535.       x = Energy(EJf[i],"K",-1);
  536.       fprintf(fp,"%f\n",x);
  537.    }
  538.    fclose(fp);
  539. }
  540.  
  541. load_data()
  542. {
  543. int i,n;
  544. char *s;
  545.  
  546.    s = (char *) malloc(100);
  547.    printf("loading from >%s< ?",fnam); fflush(stdout);
  548.    fgets(s,80,stdin);
  549.    if(s[0] != 10) {s[strlen(s)-1] = 0; strcpy(fnam,s);}
  550.    free(s);
  551.    read_data();
  552. }
  553.  
  554. read_data()
  555. {
  556. FILE *fp;
  557. int i,n;
  558. double x;
  559.  
  560.    fp = fopen(fnam,"r");
  561.    if(fp == NULL) {
  562.       fprintf(stderr,"could not open file for read: %s\n",fnam);
  563.       printf("pre calculating sums...\n");
  564.       sum_init();
  565.       printf("ready\n");
  566.       return();
  567.    }
  568.    fscanf(fp,"%lf",&gI);
  569.    fscanf(fp,"%lf",&Nq);
  570.    fscanf(fp,"%lf",&l);
  571.    fscanf(fp,"%lf",&Nel);
  572.    fscanf(fp,"%lf",&J_grund);
  573.    fscanf(fp,"%lf",&R_mittel);
  574.    fscanf(fp,"%lf",&B_core);
  575.    fscanf(fp,"%d",&NTerm);
  576.    for(i = 0; i < MAXTERM; i++) {
  577.       fscanf(fp,"%lf",&Jges[i]);
  578.       fscanf(fp,"%lf",&Lges[i]);
  579.       fscanf(fp,"%lf",&Sges[i]);
  580.       fscanf(fp,"%lf",&x);
  581.       EJf[i] = Energy(x,"K",1);
  582.    }
  583.    fclose(fp);
  584.    printf("pre calculating sums...\n");
  585.    sum_init();
  586.    printf("ready\n");
  587. }
  588.  
  589. interactive()
  590. {
  591. int i, flag;
  592. char *s,c;
  593.  
  594.    s = (char *) malloc(100);
  595.    flag = 1;
  596.    while(flag == 1) {
  597.       printf("<S>ave (data)   <M>odify (data)   <C>alculate (spectrum)  <P>lot (spectrum)\n");
  598.       printf("<L>oad (data)   <F>lag (1/T toggle)   <T>able (Beta(T))   <Q>uit\n");
  599.       printf("<G> (gj,...)\n");
  600.       fgets(s,80,stdin);
  601.       c = s[0];
  602.       switch(c) {
  603.          case 'S': case 's': save_data(); break; 
  604.          case 'L': case 'l': load_data(); break; 
  605.          case 'M': case 'm': ask_data(); break;
  606.          case 'C': case 'c': calculate_data(); break;
  607.          case 'P': case 'p': plot_spectrum(); break;
  608.          case 'F': case 'f': Flag_1dT = Flag_1dT ^ 1; break; 
  609.          case 'T': case 't': table(); break;
  610.          case 'Q': case 'q': flag = 0; break;
  611.          case 'G': case 'g': show_faktor(); break;
  612.          default: printf("unknown command >%s<\n",s);
  613.       }
  614.    }
  615.    free(s);      
  616. }
  617.  
  618. main(argc,argv)
  619. int argc;
  620. char *argv[];
  621. {
  622. int n,m,i,max;
  623. char *z;
  624.  
  625.    spc = (float *) calloc(_MAXSPCLEN,sizeof(float));
  626.    err = (float *) calloc(_MAXSPCLEN,sizeof(float));
  627.    tim = (float *) calloc(_MAXSPCLEN,sizeof(float));
  628.    EJf = (double *) calloc(MAXTERM,sizeof(double));
  629.    Lges = (double *) calloc(MAXTERM,sizeof(double));
  630.    Sges = (double *) calloc(MAXTERM,sizeof(double));
  631.    Jges = (double *) calloc(MAXTERM,sizeof(double));
  632.    sum1 = (double *) calloc(MAXTERM,sizeof(double));
  633.    sum2 = (double *) calloc(MAXTERM,sizeof(double));
  634.    sum3 = (double *) calloc(MAXTERM,sizeof(double));
  635.    z = (char *) malloc(80);
  636.    fnam = (char *) malloc(80);
  637.    strcpy(fnam,"BHF_dat");
  638.  
  639.    for(i = 0; i < MAXTERM; i++) EJf[i] = (double) (i * 1000);
  640.  
  641.    R_mittel = ReCuR(6.796,"a.u.",1);
  642.    EJf[1] = Energy(1600.0,"K",1);
  643.    B_core = 10.0;
  644.    calc_ELSJ();
  645.  
  646.    if(checkopt(argc,argv,"-r",z)) strcpy(fnam,z);
  647.    if(checkopt(argc,argv,"-T",z)) Flag_1dT = 0;
  648.    if(checkopt(argc,argv,"-Ts",z)) sscanf(z,"%lf",&Tbeg);
  649.    if(checkopt(argc,argv,"-Te",z)) sscanf(z,"%lf",&Tend);
  650.    if(checkopt(argc,argv,"-Ti",z)) sscanf(z,"%lf",&Tstep);
  651.  
  652.    read_data();
  653.    if(checkopt(argc,argv,"-go",z)) {
  654.       generate_spectrum(1);
  655.    } else {
  656.       help();
  657.       interactive();
  658.    }
  659.  
  660.    free(fnam); free(z);
  661.    free(sum1); free(sum2); free(sum3);
  662.    free(Jges); free(Sges); free(Lges);
  663.    free(EJf); free(tim); free(err); free(spc);
  664. }
  665.  
  666. SysCall(str)
  667. char *str;
  668. {
  669. int i,n,m;
  670.  
  671. #ifdef AMIGA
  672. struct FileHandle *input, *output;
  673.  
  674. input = (struct FileHandle *) Open("nil:",MODE_OLDFILE);
  675. output = (struct FileHandle *) Open("nil:",MODE_NEWFILE);
  676.    Execute(str,input,output);
  677. Close(input);
  678. Close(output);
  679. #endif
  680.  
  681. #ifdef UNIX
  682.    system(str);
  683. #endif
  684. }
  685.